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Abstract 

We modify the Kuramoto model for synchronization on complex networks by introducing a gauge 
term that depends on the edge betweenness centrality (BC). The gauge term introduces additional 
phase difference between two vertices from to ir as the BC on the edge between them increases 
from the minimum to the maximum in the network. When the network has a modular structure, 
the model generates the phase synchronization within each module, however, not over the entire 
system. Based on this feature, we can distinguish modules in complex networks, with relatively 
little computational time of O(NL), where N and L are the number of vertices and edges in the 
system, respectively. We also examine the synchronization of the modified Kuramoto model and 
compare it with that of the original Kuramoto model in several complex networks. 

PACS numbers: 89.75.-k, 89.65.-s 



Complex networks have drawn considerable attention from diverse disciplines such as so- 
ciology, information science, physics, biology and so on Many complex networks in real 
world contain modules within them, which form in a self-organized way to achieve the effi- 
ciency functionally or regionally Such modular systems can exhibit collective synchronized 
patterns within each module, not forming the global synchronization [2| as can be found in 
the cortex of neural network [3( or different synchronization transition behaviors depending 



on the patterns of inter-modular connections 



4]- 



In this Letter, we study the modular synchronization pattern generated from a modified 
Kuramoto equation (KE), which we call the gauge KE, 

d(j){t) N 

Here, <pi is the phase of vertex i, fij is the natural frequency of vertex i selected from the 
Gaussian distribution e~^ 2 / 2 / \phx , J is the overall coupling constant and Oy is the (i,j)-th 
component of the adjacency matrix, which is one when the vertices i and j are connected, 
and zero otherwise, rj is a control parameter. The extra phase term g(bij), we call the gauge 
term below, is defined as 
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where and 6 max are the minimum and the maximum edge betweenness centrality (BC) [51] 
or load |6(, respectively, in the system. Here, the edge BC or load is the amount of effective 
traffic passing through a given edge when every pair of vertices sends and receives a unit 
packet that travels along the shortest path between them. Then the gauge term g(b{j) is in 
the range from to 7r depending on the BC of edge. When rj = 0, the gauge KE recovers 
the standard KE [7| which becomes fully synchronized when J is sufficiently large. The KE 
with the extra phase of the form sin(0, — <ftj — c) (c = constant) was studied first in 
The effect of the extra phase is to destroy the synchronization. Intuitively, one expect that 
the BCs on intra-module links are smaller than those on inter-module. Thus, each module 
can be synchronized, while the entire system is not. Moreover, the gauge term induces an 
effective coupling that can be negative at the edges connecting different modules. Due to 
this negative coupling, the average phase of each module may have velocity different from 
each other. Using this property, the gauge KE can be used for module identification in 
complex networks. 
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The module identification in the context of synchronization has been studied 
These studies are inspired by the so-called dynamic clustering (DC) approach that individual 
oscillators have different levels of synchronization time owing to the heterogeneity of degree 
in network. Since vertices within modules are densely connected, they are synchronized 
more earlier than those between modules. Using this idea, the hierarchical structure can 
be detected by monitoring the temporal evolution of synchronization To identify the 
modules, however, the information of characteristic time at each hierarchical level is needed, 
which may be obtained from the spectrum of the Laplacian matrix of the system. Boccaletti 



et al., ICfl introduced another model, in which the coupling strength of the KE depends on 
the BC as where a(t) is negative. Thus, the coupling strengths across the module- 
connecting edges are weaker than those within module. a(t) is then tuned to detect the 
modules. In both methods, one needs to control the parameters such as time and a(t). 
However, our method based on Eq. ([1]) with 77 = 1 does not contain any control parameter, 
so that we can identify the modules without any prerequisite information. 

We begin to study the synchronization pattern generated from Eq. (CQ). Firstly, we apply 
the gauge KE to an ad hoc network |ll|] with a modular structure. The network is composed 
of N = 128 vertices and L = 1024 edges. Those vertices are grouped to four modules, each 
of which is of equal size. And edges are connected with probability pi n for pairs of nodes 
belonging to the same module whereas pairs belonging to different modules have edges with 
probability p ou t- By controlling the parameter pi n and p ou t we can obtain a fraction of inter- 
modular edges, z ont /(k) as we want, where z out is the mean degree of inter-modular edges 
and (k) = 2L/N is the mean degree. This ad hoc network has been used as a benchmark 



for module identification algorithms in previous studies 12]. 
We measure the order parameter defined as 
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where (• ■ ■ ) denotes the time and ensemble average. The order parameter is measured in 
the steady state. When 77 = 0, the order parameter saturates to 1 for large J, however, as 
rj is increased toward 1, it saturates at lower values as shown in Fig. QJa). This behavior 
indicates that the network is not synchronized globally. To check if the synchronization 
forms within each module, the local order parameter, defined as M. a = (\ J2f=i e%< ^ ] /N a \), is 
measured, where a is the module index, N a is the number of vertices within the module a 
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FIG. 1: (Color online) The order parameter denned over the entire network (a) and within a 
module (b) versus the coupling constant J for the ad hoc network in case of z out /(k) = 0.05. Data 
are for 77 = 0.0, 0.6, 0.7, 0.8, 0.9 and 1.0 from the top in (a). The same symbols are used for (b), 
but data for different rj collapse onto the single curve. 
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FIG. 2: (Color online) The time evolution of average phases of the four modules, distinguished by 
different symbols, for the ad hoc network with z out /(k) = 0.05 when 77 = 1.0 and J = 2.0. 

and the sum is over vertices within the module. We find that indeed the order parameter 
■M-mod reaches 1 for large J as shown in Fig. [T](b), indicating that the oscillators within the 
module are synchronized. We examine the average phase of each module as a function of 
time. As shown in Fig. [21 the modules are distinguishable by different average phases and 
average phase velocities. 

The stability of synchronization of the model is examined. Assuming the fully synchro- 
nized state of the form 0* = 0° + Qt, and linearizing Eq.(JT]), we get £j(t) = — JYuj (*ij£j(t) 



4 



-0.4 

0.2 0.4 0.6 0.8 1 
il 

FIG. 3: (Color online) The first 4 eigenvalues, Ai = 0, A2, A3 and A4, of Gij versus the parameter 
77 for the ad hoc network in case of z ont /{k) = 0.05 and J = 2.0. Data beyond rj c ~ 0.59 depend 
sensitively on time t where LOij is obtained. 

where = - 0*, = (Y,k a ^ik)Sij - a^uty and u tj = cos(0° - 4>] - Vd{bij))- 
Ai = is the trivial eigenvalue of G and the sign of other eigenvalues determines the sta- 
bility of the fully synchronized state. Due to the negative element of the coupling matrix 
G, its eigenvalues can be negative, and then the Lyapunov exponent in the linear stability 
analysis can be as well. In that case, the synchronization is no longer stable. We obtain 
Uij from cos(0j(t) — 4>j(t) — Tjgfyj)) at an arbitrary but sufficiently large t and trace out 
the eigenvalues for the ad hoc network having z out / (k) = 0.05 and plot the first 3 non-zero 
eigenvalues versus 77 in Fig. [31 A2 is positive at 77 = and decreases to zero as 77 increases 
from to r] c ~ 0.59. And increasing rj further above rj c drives the system to unstable state. 
For < 77 < r] c , the order parameter A^tot is almost 1 in the steady state, whereas A^tot has 
a smaller constant value for rj > r\ c . In many cases, they actually oscillates in time before the 
time average due to disparate group velocities of the modules as shown in Fig. [2j The curve 
fitting of A2 in the vicinity of 77 = r/ c shows A2 oc (r? c — 77) l l 2 . The square-root singularity of 
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A2 near the stability edge is the signature of the saddle-node bifurcation 

We introduce how to identify modules with the gauge KE. To this end, we take the 
following steps: 

i) We apply the gauge KE ([T]) to all oscillators with a sufficiently large coupling constant 
J. The phases {4>i(t)} of each oscillator are obtained in the steady state. 

ii) We measure the phase similarity defined as = ([1 + cos(0»(t) — <pj(t))]/2) for each 



connected pair of oscillators The brackets are the average over different times, 

natural frequencies {fli}, and initial random phases {0,(0)}. 

Hi) From the empty state, where all edges are absent, we add edges (i, j) one by one that 
are chosen following the descending order of CV,-. 

Clusters after the step Hi) are regarded as modules. The edges that existed originally, 
but not connected yet until the step Hi) are regarded as inter-modular edges. 

iv) We repeat the step Hi) until the modularity of the system becomes maximum. The 
modularity Q is defined as 

Q = ^e aa - a 2 a , (4) 

a 

where a a = e a/ 3, and e a p is the fraction of edges that connect the vertices belonging 



to the modules a and (3 [H 



To test the performance of our algorithm, we measure the mutual information on several 
networks, defined as 

[ ' } Yf=i Ni log(f ) + M log(^) (5j 

where M = 4 is the number of preassigned modules and M' is the number of detected 
modules. N- is the number of vertices belonging to the i-th preassigned and the j-th 



detected modules, N { = J2j N l and N 3 = J2i N l [12]- 

Fig. H] shows the mutual information measured on the ad hoc network as a function of 
z ut/(k) for several module detecting algorithms. The performance of our algorithm is not 
better than those of the Potts model and the simulated annealing (SA) Q, [l5]. Even though 
they are better in performance, if we count for their long computation time, then ours may be 



useful practically. The performance of opinion changing rate model (OCR) algorithms [10| is 
somewhat better, however, it requires an extra task of parameter tuning, so that ours is easier 
to implement. Since our algorithm shares with the Girvan-Newman (GN) algorithm [16] the 
idea of clustering based on BC, the performances of the two algorithms are close to each 
other. However, since ours calculates the BC on each edge only once, whereas the GN 
algorithm does it repeatedly for each disconnected cluster, the computational time can be 
reduced drastically from 0(NL 2 ) to O(NL). The performance of our algorithm is better 
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FIG. 4: (Color online) The mutual information versus z on t/{k), the fraction of inter-modular edges 
per mean degree for the ad hoc network. See the text for abbreviations. 
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FIG. 5: (Color online) The dendrogram based on the phase similarity between connected pairs of 
vertices for the hierarchical network with three levels. 



than that of the Clauset-Newman-Moore (CNM) algorithm |l7j . which runs in C(iVln 2 N) 
for sparse graphs. 

Secondly, we apply our algorithm to the hierarchical network proposed by Ravasz and 



Barabasi 18|. When the number of levels is two, the modules are well selected in a similar 
way as in Fig. 3 of Ref. For the three level case, the dendrogram constructed by our 
method is shown in Fig. [5j Here, the hub at the second level is grouped with one of the four 
identical modules connected to it in the second level. 

Thirdly, we apply the gauge KE to Erdos-Renyi (ER) random networks and scale-free 
(SF) networks with no modular structures to see the network structure dependences. The SF 
network is generated using the static model 6J. The order parameter ([3]) behaves differently 
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FIG. 6: (Color online) The order parameter versus the coupling constant J for the ER (a) and the 
SF network with the degree exponent 3.5 (b). The data are for the cases of rj = 0,0.2,0.4,0.6,0.8 
and 1.0 from the top. 
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FIG. 7: (Color online) The phase difference across the edge with the maximum BC. 

for the two networks. For the ER network, the saturated value of the order parameter 
decreases from 1 to as 77 increases from to l(Fig. However, for the SF network, 

the order parameter does not decrease to 0, but w 0.7 even if 77 reaches l(Fig. E(b)). To 
study the origin of the different behaviors, we measure the phase difference A<p across the 
edge with the maximum BC. In most cases, one end of the edge is the hub. For the ER 
network, its change with time is large running from — 7r to 7r as shown in Fig. For the SF 
network, it stays around a smaller value in short intervals. Such difference is rooted from 
the following. For the SF network, the hub has large degree, so that the probability to form 
a triangle including the miximum BC edge is larger for the SF network than for the ER 
network, provided that the mean degree of the system is the same. Owing to such short 
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loops, the phase difference across the maximum BC edge is small for the SF network, and 
large for the ER network. The overall order parameter is close to zero for the ER network. 

In summary, we have introduced a gauge KE in which the gauge term depends on the 
edge BC. The gauge term drives the phase difference between the two vertices of an edge 
from to 7r as the BC across the edge increases. As a result, the phase difference of two 
oscillators belonging to different modules is large, however, it is small across the edges within 
modules. Thus, the model generates the phase synchronization within each module, however, 
it does not globally. Measuring the phase similarity between two connected oscillators, we 
constructed the dendrogram and identified the modules. Such module detecting method 
works efficiently 
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